Integral equation formulations are methodologies that consists on transforming the differential operators into integral operators exactly incorporating the boundary conditions. This work extends the use of integral equations formulations (IEF) for the development of numerical schemes for generalized reaction-diffusion systems. Due to the property of exact incorporation of the boundary conditions, the numerical schemes resulting do not require approximations at the boundaries, which leads to numerical schemes with global orders approximation O(h²). IEF schemes are characterized by their systematic methodology and its their mathematical structure that allows easy physical interpretation of the processes involved in the reaction-diffusion systems. To evaluate the ability of numerical approximation of the proposed scheme different conditions that arise in the study of catalytic pellets were considered. The numerical simulations carried out show that the IEF exhibit better numerical approximations than classic finite differences schemes (FD).